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The inelastic Boltzmann equation for a granular gas is applied to spatially inhomogeneous states 
close to the uniform shear flow. A normal solution is obtained via a Chapman-Enskog-like expan- 
sion around a local shear flow distribution. The heat and momentum fluxes are determined to 
first order in the deviations of the hydrodynamic field gradients from their values in the reference 
state. The corresponding transport coefficients are determined from a set of coupled linear integral 
^ ' equations which are approximately solved by using a kinetic model of the Boltzmann equation. The 

main new ingredient in this expansion is that the reference state / ( - - ) (zeroth-order approximation) 
' retains all the hydrodynamic orders in the shear rate. In addition, since the collisional cooling 

cannot be compensated locally for viscous heating, the distribution /^°^ depends on time through 
0^ , its dependence on temperature. This means that in general, for a given degree of inelasticity, the 

complete nonlinear dependence of the transport coefficients on the shear rate requires the analysis 
of the unsteady hydrodynamic behavior. To simplify the analysis, the steady state conditions have 
been considered here in order to perform a linear stability analysis of the hydrodynamic equations 
with respect to the uniform shear flow state. Conditions for instabilities at long wavelengths are 
identified and discussed. 
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I. INTRODUCTION 



The understanding of granular systems still remains a topic of interest and controversy. Under rapid flow conditions, 
they can be modeled as a fluid of hard spheres dissipating part of their kinetic energy during collisions. In the 
fyT) • simplest model, the grains are taken to be smooth so that the inelasticity of collisions is characterized through a 
constant coefficient of normal restitution a < 1. Energy dissipation has profound consequences on the behavior of 
these systems since they exhibit a rich phenomenology with many qualitative differences with respect to molecular 
systems. In particular, the absence of energy conservation yields subtle modifications of the conventional Navier- 
Stokes equations for states with small gradients of the hydrodynamic fields. The dependence of the corresponding 
transport coefficients on dissipation may be determined from the Boltzmann kinetic equation conveniently modified 
to account for inelastic binary collisions 0, Q • The idea is to extend the Chapman-Enskog method Q to the inelastic 
case by expanding the velocity distribution function around the local version of the homogeneous cooling state, 
namely, a homogeneous state whose dependence on time occurs only through the temperature. In the first order of 
the expansion, explicit expressions for the transport coefficients as functions of the coefficient of restitution have been 
obtained in the case of a single gas Q as well as for granular mixtures showing good agreement the analytical 
results with those obtained from Monte Carlo simulations @. 

Although the Chapman-Enskog method can be in principle applied to get higher orders in the gradients (Burnett 
and super-Burnett corrections,. . .), it is extremely difficult to evaluate those terms especially for inelastic systems. In 
addition, questions about its convergence remains still open 0. This gives rise to the search for alternative approaches 
to characterize transport for strongly inhomogeneous situations (i.e., beyond the Navier-Stokes limit). One possibility 
is to expand in small gradients around a more relevant reference state than the (local) homogenous cooling state. For 
example, consider states near a shearing reference steady state such as the so-called uniform (simple) shear flow (USF) 
0. Such an application of the Chapman-Enskog method to a nonequilibrium state requires some care as recently 
discussed in Ref. 0- The USF state is probably the simplest flow problem since the only nonzero hydrodynamic 
gradient is du x /dy = a = const, where u is the flow velocity and a is the constant shear rate. Due to its simplicity, 
this state has been widely used in the past both for elastic and inelastic gases to shed light on the complexities 
associated with the nonlinear response of the system to the action of strong shearing. However, the nature of this 
state for granular systems is different from that of the elastic fluids since the source of energy due to the macroscopic 
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imposed shear field drives the granular system into rapid flow and a steady state is achieved when the amount of 
energy supplied by shearing work is balanced by the lost one due to the inelastic collisions between the particles. As a 
consequence, in the steady state the reduced shear rate a* cx a/y/T (which is the relevant nonequilibrium parameter 
of the problem) is not an independent quantity but becomes a function of the coefficient of restitution a. This means 
that the quasielastic limit (a — * 1) naturally implies the limit of small shear rates (a* « 1) and vice versa. The 
study of the rheological properties of the USF state has received a great deal of attention in recent years in the case 
of monocomponent [IS[illil[ll[ll[ll[ll[Ii[ll[ll|23and multicomponent systems [H |H M, |H |H E3 . 

The aim of this paper is to determine the heat and momentum fluxes of a gas of inelastic hard spheres under 
simple shear flow in the framework of the Boltzmann equation. The physical situation is such that the gas is in 
a state that deviates from the simple shear flow by small spatial gradients. The starting point of this study is a 
recent approximate solution of the Boltzmann equation which is based on Grad's method (2^, |U H3 • In spite of 
this approach, the relevant transport prop erties obtained from this solution compare quite well with Monte Carlo 
simulations even for strong dissipation |l8tl22ll28| ]. showing again the reliability of Grad's approximation to compute 
the lowest velocity moments of the velocity distribution function. Since the system is slightly perturbed form the 
USF, the Boltzmann equation is solved by applying the Chapman-Enskog method around the (local) shear flow state 
rather than the (local) homogeneous cooling state. This is the main feature of this expansion since the reference state 
is not restricted to small values of the shear rate. One important point is that, for general small deviations from the 
shear flow state, the zeroth-order distribution is not a stationary distribution since the collisional cooling cannot be 
compensated locally for viscous heating. This fact gives rise to new conceptual and practical difficulties not present 
in the previous analysis made for elastic gases to describe transport in thermostatted shear flow states |2^. Due to 
the difficulties involved in this expansion, here general results will be restricted to particular perturbations for which 
steady state conditions apply. In the first order of the expansion, the generalized transport coefficients are given 
in terms of the solutions of linear integral equations. To get explicit expressions for these coefficients, one needs to 
know the fourth-degree moments of USF. This requires to consider higher-order terms in Grad's approximation for 
the reference distribution function, which is quite an intricate problem. In order to overcome such difficulty, here I 
have used a convenient kinetic model |3lj that preserves the essential properties of the inelastic Boltzmann equation 
but admits more practical analysis. The mathematical and physical basis for this model as a good representation of 
the Boltzmann equation has been discussed in Ref. |3lj . In particular, it is worth noting that the results derived from 
this model coincides with those given from the Boltzmann equation at the level of the rheological properties 0, . 
Furthermore, recent computer simulation results |28( have also shown good agreement between the kinetic model and 
the Boltzmann equation for the fourth-degree moments, covering this agreement a wide range of values of dissipation 
(say, for instance, a > 0.5). This good agreement extends that previously demonstrated for Couette flow in dilute 
gases [3(J and for USF in dense systems [2(| and shows the reliability of the kinetic model to capture the main trends 
of the Boltzmann equation, especially those related to transport properties. 

The knowledge of the above generalized transport coefficients allows one to determine the hydrodynamic modes 
from the associated linearized hydrodynamic equations. This is quite an interesting problem widely analyzed in the 
literature. As noted by the different molecular dynamics experiments carried out for the USF problem |TaL l25t f3^| . it 
becomes apparent the development of inhomogeneities and formation of clusters as the flow progresses. Consequently, 
the USF state is unstable for long enough wavelength spatial perturbations. In order to understand this phenomenon, 
several stability analysis have been undertaken [23, Hj, [2H, |3(| ■ Most of them are based on the Navier-Stokes 
equations |.3.1 l34i |35j and, therefore, they are limited to small velocity gradients, which for the USF problem means 
small dissipation. Another alternative has been to solve the Boltzmann equation by means of an expansion in a set of 
basis functions |36l l37l ] . The coefficients of this expansion are then determined by using also an expansion in powers 
of the parameter e = \/l — a 2 , which is assumed to be small. All these analytical results have shown that the USF 
becomes unstable for certain kind of disturbances. My approach is different from previous works since the conditions 
for stability arc obtained from a linear stability analysis involving the transport coefficients of the perturbed USF state 
instead of the usual Navier-Stokes coefficients. Furthermore, the analysis is not restricted to the low-dissipation limit 
since the reference state goes beyond this range of values of a. Two different perturbations to the reference state have 
been considered here: (i) perturbations along the velocity gradient (y direction) only and (ii) perturbations along 
the vorticity direction (z direction) only. The results show that the USF is linearly stable in the first case while it 
becomes unstable in the second case. These results agree qualitatively with those previously derived |33l l35l | in the 
context of the Navier-Stokes description. On the other hand, at a quantitative level, the comparison carried out here 
shows significant differences between the Navier-Stokes description and the present results as the collisions become 
more inelastic. In addition, our results also confirm that the instability is confined to long wavelengths (small wave 
numbers) and so it can be avoided for small enough systems. 

The plan of the paper is as follows. In Sec. ^ the Boltzmann kinetic equation is introduced and a brief summary 
of relevant results concerning the USF problem is given. In Sec. lIIII the problem we arc interested in is described and 
the set of generalized transport coefficients characterizing the transport around USF is defined. Explicit expressions 
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for these coefficients are provided in Sec. IIVI bv using a kinetic model of the Boltzmann equation. The details of the 
calculations are displayed along several Appendices. Section [V] is devoted to the linear stability analysis around the 
steady USF state and presents the form of the hydrodynamic modes. The paper is closed in Sec. I VII with a discussion 
of the results obtained here. 

II. BOLTZMANN KINETIC EQUATION AND UNIFORM SHEAR FLOW 

Let us consider a granular gas composed by smooth spheres of mass m and diameter a. The inelasticity of collisions 
among all pairs is accounted for by a constant coefficient of restitution < a < 1 that only affects the translational 
degrees of freedom of grains. In a kinetic theory description all the relevant information on the state of the system 
is given by the one-particle velocity distribution function /(r, v, t). At low density the inelastic Boltzmann equation 
HHI gives the time evolution of /(r, v, t). In the absence of an external force, it has the form 

I- +v.vWv,t)=J[v|/(t), /(*)], (1) 



where the Boltzmann collision operator is 

J[vi\f,f] = <r 2 J dv 2 j dae(*-g)(*-g) 

x [ a - 2 /(r,vi)/(r,v^i)-/(r, Vl ,i)/(r,v 2 ,<)] . (2) 

Here, <r is a unit vector along their line of centers, is the Heaviside step function, and g = vi — v 2 is the relative 
velocity. The primes on the velocities denote the initial values {v' 1; v 2 } that lead to {vi,v 2 } following a binary 
collision: 

vl = vj - - {l + a- 1 ) (ff-g)S, V 2 = V 2 + - (l + a" 1 ) • g)£ (3) 
The first five velocity moments of / define the number density 

n(r,t) = J dv/(r,v,t), (4) 



the flow velocity 



and the granular temperature 



U(r ' <) = ^)/ dVV/(r ' V ' i} ' 

T(r,t) = — ^— f dv^ 2 (r,t)/(r,v,t), 
3ra(r,i) J 



(5) 



(6) 



where V(r, f) = v — u(r, t) is the peculiar velocity. The macroscopic balance equations for density n, momentum mu, 
and energy |nT follow directly from Eq. by multiplying with 1, mv, and ^mv 2 and integrating over v: 

D t n + nV • u = , (7) 



D t Ui + {mn)- l VjPij = , (8) 

D t T+-^(V -q + PijVjiii) = -CT, (9) 

where Dt = dt + u • V. The microscopic expressions for the pressure tensor P, the heat flux q, and the cooling rate ( 
are given, respectively, by 



P(r,t)= / dvmWf(r,v,t), (10) 
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q(r,t) = J dvimy 2 V/(r,v,i), (11) 

«'•*) = "3^Tm/ dV ^ 2j[r ' V|/W] ' (12) 

We assume that the gas is under uniform (or simple) shear flow (USF). This idealized macroscopic state is char- 
acterized by a constant density, a uniform temperature and a simple shear with the local velocity field given by 

Ui = a^rj, a l j=a5 lx 5 jy , (13) 

where a is the constant shear rate. This linear velocity profile assumes no boundary layer near the walls and is 
generated by the Lee-Edwards boundary conditions |38[ . which are simply periodic boundary conditions in the local 
Lagrangian frame moving with the flow velocity. For elastic gases, the temperature grows in time due to viscous 
heating and so a steady state is not possible unless an external (artificial) force is introduced Q. However, for 
inelastic gases, the temperature changes in time due to the competition between two (opposite) mechanisms: on the 
one hand, viscous (shear) heating and, on the other hand, energy dissipation in collisions. A steady state is achieved 
when both mechanisms cancel each other and the fluid autonomously seeks the temperature at which the above 
balance occurs. Under these conditions, in the steady state the balance equation © becomes 



aP*v = -oCP, (14) 



where p — nT is the hydrostatic pressure. Note that for given values of the shear rate a and the coefficient of 
restitution a, the relation (|14|) gives the temperature T in the steady state as a unique function of the density n. 

The USF problem is perhaps the nonequilibrium state most widely studied in the past few years both for granular 
and conventional gases At a microscopic level, it becomes spatially homogeneous when the velocities of the 

particles are referred to the Lagrangian frame of reference co-moving with the flow velocity u [3|j. Therefore, the 
one-particle distribution function adopts the uniform form, /(r, v) — * /(V), and the Boltzmann equation reads 

-aV y ^-J(V) = J[V\fJ]. (15) 

This equation is invariant under the transformations 

V z ^-V z , (V x ,V y ) -> -(V x ,V y ), (V x ,a)^(-V x ,-a). (16) 

The elements of the pressure tensor provide information on the relevant transport properties of the USF problem. 
These elements can be obtained by multiplying the Boltzmann equation l|15|) by mViVj and integrating over V. The 
result is 

auPjt + a jt P u = m J dVV t V r J[V\fJ] 

= M r (17) 

The exact expression of the collision integral Ajj is not known, even in the elastic case. However, a good estimate can 
be expected by using Grad's approximation: 



/(V) -> / (V) 

where /o(V) 



(18) 



/ (V) = n(m/27rT) 3/2 exp(-mV r2 /27 1 ) (19) 

is the local equilibrium distribution function. When Eq. I|18(l is substituted into the definition of Ay and nonlinear 
terms in Pij/nT — 5^ are neglected, one gets 

A y - = -v \fi (R tJ - pS i:j ) + CPA , (20) 
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where 



is an effective collision frequency, 
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C* = ^ = ^(l-« 2 ), (22) 



is the dimensionless cooling rate evaluated in the local equilibrium approximation and 

The set of coupled equations for P,j can be now easily solved when one takes into account the approach (|2L)fl . The 
expressions for the reduced elements P*j = -P^/p are 

P* = 3 — 2 P* P* = P* = - P* = - n* (24) 

« vv yy zz p _|_ ^* ' ^ _|_ £*)2 ' \ ' 

where the (reduced) shear rate a* = a/i> is given by 

/ Q /•* 

(25) 




The expression (|25ll clearly indicates the intrinsic connection between the (reduced) velocity gradient and dissipation 
in the system. In fact, in the elastic limit (a = 1, which implies a* =0), the equilibrium results of the ordinary gas 
are recovered, i.e., Ph = <5y. This means that a (or a*) can be considered as the relevant nonequilibrium parameter of 
the problem. The analytical results given by Eqs. I)24|) and H25f) agree quite well H^,[2]j with Monte Carlo simulations 
of the Boltzmann equation [22L l28l| , even for strong dissipation. 



III. SMALL PERTURBATIONS FROM THE UNIFORM SHEAR FLOW: TRANSPORT COEFFICIENTS 

In general, the USF state can be disturbed by small spatial perturbations. The response of the system to these 
perturbations gives rise to additional contributions to the momentum and heat fluxes, which can be characterized by 
generalized transport coefficients. This section is devoted to the study of such small perturbations. 

In order to analyze this problem we have to start from the Boltzmann equation with a general time and space 
dependence. Let Uo = a • r be the flow velocity of the undisturbed USF state. Here, the only nonzero element 
of the tensor a is a,j = a8i X 8j y . In the disturbed state, however the true velocity u is in general different from 
u since u = u + <5u, Su being a small perturbation to u . As a consequence, the true peculiar velocity is now 
c = v u = V Su, where V = v — u . In the Lagrangian frame moving with u , the Boltzmann equation can be 
written as 



t T/ 9 

Ot oV, 



f - aV yJ^rf + ( V + u o) ' V/ = J[V|/, /], (26) 



where here the derivative V/ is taken at constant V. The corresponding macroscopic balance equations associated 
with this disturbed USF state follows from the general equations © when one takes into account that u = Uo + <5u. 
The result is 

d t n + u • Vn = -V • (nSu), (27) 
d t Su + a • 5u + (u + Su) ■ VSu = -(mn)- 1 ^ ■ P, (28) 



^nd t T + ^n(u + <5u) ■ VT + aP xy + V ■ q + P : V<5u = ~pC, (29) 



6 



where the pressure tensor P, the heat flux q and the cooling rate £ are defined by Eqs. (|10(l - (|12|l . respectively, with 
the replacement V — > c. 

We assume now that the deviations from the USF state are small, which means that the spatial gradients of the 
hydrodynamic fields 

A(r,t) = {n(r,t),T(r,t),*u(r,t)} (30) 

are small. Under these conditions, a solution to the Boltzmann equation (|26|l can be obtained by means of a general- 
ization of the conventional Chapman-Enskog method |3j where the velocity distribution function is expanded about a 
local shear flow reference state in terms of the small spatial gradients of the hydrodynamic fields relative to those of 
USF. This type of Chapman-Enskog-like expansion has been considered in the case of elastic gases to get the set of 
shear-rate dependent transport coefficients 0, in a thermostatted shear flow problem and it has also been recently 
considered ||| in the context of inelastic gases. 

To construct the Chapman-Enskog expansion let us look for a normal solution of the form 

f(r,y,t)=f(A(r,t),V). (31) 

This special solution expresses the fact that the space dependence of the reference shear flow is completely absorbed 
in the relative velocity V and all other space and time dependence occurs entirely through a functional dependence 
on the fields A(r,t). The functional dependence can be made local by an expansion of the distribution function in 
powers of the hydrodynamic gradients: 

/(r, V, t) = (A(r, t), V) + (A(r, t), V) + • • • , (32) 

where the reference zeroth-order distribution function corresponds to the USF distribution function but taking into 
account the local dependence of the density and temperature and the change V — * V — Su(r, t) [see Eqs. (|D3|I and 
(|D4|I for the explicit form of /^°^ in the steady state given by a kinetic model of the Boltzmann equation]. The 
successive approximations /' fe ' are of order k in the gradients of n, T, and Su but retain all the orders in the shear 
rate a. This is the main feature of this expansion. In this paper, only the first order approximation will be considered. 
More details on this Chapman-Enskog-like type of expansion can be found in Ref . y| ■ 

The expansion (|32|l yields the corresponding expansion for the fluxes and the cooling rate when one substitutes 
into their definitions (|TUj| - [jT2T) : 

P=P (0) + p(i) + ... t q=q (0) +q (i) + ... > ( = (W +(W + ■■■ . (33) 

Finally, as in the usual Chapman-Enskog method, the time derivative is also expanded as 

d t = 9 t (0) + d ( t 1] + d[ 2) + • • • , (34) 

where the action of each operator is obtained from the hydrodynamic equations (|27|) - (|29|l . These results provide 
the basis for generating the Chapman-Enskog solution to the inelastic Boltzmann equation l|2tj|) . 

A. Zeroth-order approximation 

Substituting the expansions IK>2[I and l|34|l into Eq. I|26f) . the kinetic equation for f(°' is given by 

dl 0) f m -aV y ^J^=J[V\f(°\f(% (35) 
To lowest order in the expansion the conservation laws give 



#n = 0, #T=-A aP (0)_ rC (0) ) (36) 
3n 

d^Sui+aijSuj = 0. (37) 

As said before, for given values of a and a, the steady state condition (|14|) establishes a mapping between the density 
and temperature so that every density corresponds to one and only one temperature. Since the density n(r,t) and 
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temperature T(r, t) are specified separately in the local USF state, the viscous heating only partially compensates for 
the collisional cooling and so, d^T ^ 0. Consequently, the zeroth-order distribution /(°) depends on time through 
its dependence on the temperature. Because of the steady state condition 114(1 does not apply in general locally, the 
reduced shear rate a* — a/v(n,T) depends on space and time so that, a* and a must be considered as independent 
parameters for general infinitesimal perturbations around the USF state. Since is a normal solution, then 



a m fW> 



On 



9/ (0) a( o) 
dT ^ OSu. ' 



d">6 Ui 



= -(|- a p(o) + rc(°)U^ + ^4^, 



(38) 



where in the last step we have taken into account that / ( °) depends on Su only through the peculiar velocity c. 
Substituting ((38(1 into l(35|) yields the following kinetic equation for /(°); 



AaPW +TC (0) ) ^/ (0) ~ «St|-/ (0) = ^[V|/ (0) ,/ ( °]. 



(39) 



The zeroth-order solution leads to q(°) — 0. On the other hand, to solve Eq. I(39|l one needs to know the temperature 
dependence of the zeroth momentum flux Pxy . A closed set of equations for P*- -* is obtained when one considers 
Grad's approximation ((181) : 
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p, 



(0) 



Mi. 



(0) 



(40) 



where 



c* 



£(0) 
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(41) 



The steady state solution of Eq. 1(40(1 is given by Eqs. 1(24(1 and ((25(1 . However, in general the equations 1(40(1 must 
be solved numerically to get the dependence of the zeroth-order pressure tensor P^ ' (T) on temperature. A detailed 
study on the unsteady hydrodynamic solution of Eqs. 1(401) has been carried out in Ref. [^J. In what follows, P^{T) 



will be considered as a known function of T. 



B. First-order approximation 



The analysis to first order in the gradients is worked out in Appendix^ Only the final results are presented in 
this Section. The distribution function /W is of the form 



= X„ • Vn + X T • VT + X u : VSu, 



(42) 



where the vectors X„ and and the tensor X„ are functions of the true peculiar velocity c. They are the solutions 
of the following linear integral equations: 



3n 



d 



X,, 



T 



2n 



(1 - nd n )P£> - C 



(«) 



X T A = Yn 



(43) 



AftPi") + TC (0) ) d T + fTid T P^) + § C <°> + ac y £- - C 



dc x 



Xta = Yt', 



(44) 



X„ 



aSk v X UiX £ — (u.keTdrf^ — Y u 



(45) 
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where Y n (c), Y T (c), and Y„(c) are denned by Eqs. HA9p - ljAll|l . respectively, and ( u ^ k e is denned by Eq. l|A14jl , While 
the Y functions are given in terms of the reference state distribution f^°\ C, u ,ki is a functional of the unknown X u ^g- 
In addition, £ is the linearized Boltzmann collision operator around the reference state 



LX = 



(j[fW,X] + J[X, /(»)]). (46) 



A good estimate of Cu,ki can be obtained by expanding X u ^e in a complete set of polynomials (for instance, Soninc 
polynomials) and then truncating the series after the first few terms. In practice, the leading term in these expansions 
provides a very accurate result over a wide range of dissipation. This contribution has been obtained in Appendix IbI 
and is given by Eq. (|B9(1 . 

With the distribution function /W determined by (|42|l . the first-order corrections to the fluxes are 

4" - -*m& («) 

(i) dT dn . 

« "■'•^ 11 -7^: (48) 

where 

Vijkt = ~ J dcmCiCjX uM (c), (49) 

Kij = - J dc^-c 2 CiX T j(c), (50) 

J dcyC 2 Ci X nJ (c). (51) 

Upon writing Eqs. (|47|l - (|51|l use has been made of the symmetry properties of X n j Xx,i and X Ut ij. In general, the 
set of generalized transport coefficients rjijkt-, Kij, and /iy are nonlinear functions of the coefficient of restitution a 
and the reduced shear rate a*. The anisotropy induced in the system by the shear flow gives rise to new transport 
coefficients, reflecting broken symmetry. The momentum flux is expressed in terms of a viscosity tensor rjijki{a) of 

rank 4 which is symmetric and traceless in ij due to the properties of the pressure tensor Pa- . The heat flux is 
expressed in terms of a thermal conductivity tensor Kij (a) and a new tensor /i.y(a). 

C. Steady state conditions 

As shown in the above subsections, the evaluation of the complete nonlinear dependence of the generalized transport 
coefficients on the shear rate and dissipation requires the analysis of the hydrodynamic behavior of the unsteady 
reference state. This involves the corresponding numerical integrations of the differential equations obeying the 
velocity moments of the zeroth-order solution. This is a quite intricate and long problem. However, given that here 
we are mainly interested in performing a linear stability analysis of the hydrodynamic equations with respect to the 
steady state, we want to evaluate the transport coefficients in this special case. As a consequence, d[ 0) T = and so 



Hij 



the condition 

applies. In Eq. j52jl. it is understood that a* and P* y = P^y /p are evaluated in the steady state, namely, they are 
given by Eqs. I|24|l and (|25|l . respectively. A consequence of Eq. (|52|) is that the first term on the left hand side of 
the integral equations 143ll -145 [) vanishes. In addition, the dependence of the pressure tensor P^- on density and 



temperature occurs explicitly through p — nT and through its dependence on a* . In this case, the derivatives d n pjj 
and drPjj can be written more explicitly as 



ndnPW = nd nP P* 3 (a*)=p\l-a*-^\ P* (a*), (53) 
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Td T P§ ] = TdrpP^a*) =p(l- \a*^j P^a*). 



(54) 



The dependence of P*j on a* near the steady state is determined in the Appendix IU1 so that, all the terms appearing 
in the integral equations are explicitly known in the steady state. Under the above conditions, Eqs. l(15 j) -l|15 )l become 

-acy-^- + C^j X n>i + y^( P *v + a * d *aKy)XT,i = Y n<i , (55) 
(-acy A - (P* y - a*d a ,P* y ) + C^j X T>i = Y T4 , (56) 

-acy-r^- + CJ X u .kt - a8kyX u . x i - Cu,ki T 9rf^ = Y Ut n, (57) 

where it is understood again that in Eqs. H55(l - (|57|l all the quantities are evaluated in the steady state. Henceforth, I 
will restrict my calculations to this particular case. 

Given that in the steady state the coefficient of restitution and the reduced shear rate are coupled, the usual 
Navier-Stokes transport coefficients for ordinary gases are recovered for elastic collisions (a* = 0). Thus, when a — > 1 
the coefficients become 



(58) 



where 770 = pjv and kq = 15?7o/4m are the shear viscosity and thermal conductivity coefficients given by the (elastic) 
Boltzmann equation. 

IV. RESULTS FROM A SIMPLE KINETIC MODEL 

The explicit form of the generalized transport coefficients /iy , and rjij ^ requires to solve the integral equations 
(|55H -f l|57|l . Apart from the mathematical difficulties embodied in the Boltzmann collision operator £, the fourth- 
degree velocity moments of the distribution /(°) are also needed to determine fiij and Kij and they are not provided 
in principle by the Grad approximation. Nevertheless, an accurate estimate of these moments from the Boltzmann 
equation is a formidable task since it would require at least to include the fourth-degree moments in Grad's solution. 
In this case, to overcome such difficulties it is useful to consider a model kinetic equation of the Boltzmann equation. 
As for elastic collisions, the idea is to replace the true Boltzmann collision operator with a simpler, more tractable 
operator that retains the most relevant physical properties of the Boltzmann operator. Here, I consider a kinetic 
model based on the well-known Bhatnagar-Gross-Krook (BGK) Q for ordinary gases where the operator J[f, f] 
is 

J[f, /] - -/M/ -/<>) + !!-• ( c /) • ( 59 ) 

Here, v and (3 are given by Eqs. (|21|l and (|23|l . respectively, /o is the local equilibrium distribution (|19|l and £ is the 
cooling rate defined by Eq. (|12|l . As said before, an estimate of £ to first order in the gradients has been derived 
in Appendix [B] In general, the quantity (3 can be considered as an adjustable parameter to optimize the agreement 
with the Boltzmann equation. In this paper, 8 has been chosen to reproduce the true Navier-Stokes shear viscosity 
coefficient of an inelastic gas of hard spheres [j]. A slightly different choice for 8, namely = (1 + a) /2, is considered 
in Ref. 0. 

By taking moments with respect to 1, c and c 2 , the model kinetic equation 159|) yields the same form of the 
macroscopic balance equations for mass, momentum, and energy, Eqs. Q - ©? as those given from the Boltzmann 
equation. When a = 1, then 8 = 1,(^0 and so the kinetic model (|59|l reduces to the BGK equation whose utility 
to address complex states not accessible via the Boltzmann equation is well-established for elastic gases 0- In the 
case of granular gases, it is easy to show that the kinetic model leads to the same results for the pressure tensor in the 
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FIG. 1: Fourth-degree velocity moment (c 4 ) relative to its local equilibrium value as a function of the coefficient of restitution. 
Solid line is the prediction of the kinetic model while the symbols are simulation results |2S| . 




0.5 0.6 0.7 0.8 0.9 1.0 



a 

FIG. 2: Plot of the reduced coefficients (a) y,y V and (6) \i* xy as a function of the coefficient of restitution a. 



USF problem as those given from Grad's solution to the Boltzmann equation, Eqs. Ij24(l (|25|l . This result, along with 
those of Refs. [21| and [30j . confirms the reliability of the kinetic model for granular media as well. A summary of the 
USF results derived from the kinetic model is provided in Appendix iDl In particular, beyond rheological properties, 
recent computer simulations |28( have confirmed the accuracy of the kinetic model to capture the dependence of the 
fourth-degree velocity moments (whose expressions are needed to get the coefficients /Uy and K%j on dissipation in the 
USF state. To illustrate it, in Fig. ^we plot the fourth-degree moment 

(c 4 ) =Jdc c\f(c) (60) 

relative to its local equilibrium value (c 4 )o = lhnT 2 /m? . The symbols refer to the numerical results obtained from the 
DSMC method |2^. It is quite apparent that the analytical results agree well with simulation data (the discrepancies 
between both results are smaller than 3%), showing again that the reliability of the kinetic model goes beyond the 
quasielastic limit. 

Let us consider the perturbed USF problem in the context of the kinetic model. By using the model l|59[). the 
integral equations (|55|I - H57|I still apply with the only replacement 

CX^^X-^-—-(cX), (61) 

in the case of X n ^ and Xx,i and 

in the case of X u i j. In the above equations, C,^ is the zeroth-order approximation to C which is given by Eq. I|41|l . 
With the changes l|61|l and i|62|) all the generalized transport coefficients can be easily evaluated from Eqs. ()55f) - (|57|l . 
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FIG. 4: Plot of the reduced coefficients (a) ri yyyy , (b) rjxyyy, (c) V*z yy , an d (d) r\% xyy as a function of the coefficient of restitution 
a. 



Details of these calculations are also given in Appendix El a more complete listing can be obtained on request from 
the author. 

The dependence of the generalized transport coefficients on the coefficient of restitution a is illustrated in Figs. [3 El 
and H for the (reduced) coefficients n* xy , K* y , r] xxyy , ri* yyy , r]* zyy , and rf xyyy . Here, ^ = ti^/Tkq, k*j = «y/«0 
and r)*j kl = rnjkt/r)a, where rjo = p/v and ko — 5r]o/2m are the elastic values of the shear viscosity and thermal 
conductivity coefficients given by the BGK kinetic model. In general, we observe that the influence of dissipation on 
the transport coefficients is quite significant. 

With all the transport coefficients known, the new constitutive equations (|4Tp and (|48|l are completed and the 
corresponding set of closed hydrodynamic equations H27|) - (|29[) can be derived. They are given by 

d t n + u ■ Vn + V ■ (nSu) = 0, (63) 



1 d /^( ) d8u k 



d t Sui + aijSuj + (u + Su) ■ S/Sui + —g^r: y p ij ~ Vijki-Q^- ) = 0, (64) 

3 3 dSu 

— ndtT H — n(uo + 5u) ■ VT — ar) xy ij — 



3 „ „ 3 „ „ dSui 



-^(--nTC^— (65) 

(2) 

Note also that consistency would require to consider the term aPxy which is of second order in gradients and so, 
it should be retained. Given that this would require to determine the second order contributions to the fluxes, this 
term will be neglected in our study. An important feature of our linearized hydrodynamic equations is that they are 
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not restricted to small values of the (reduced) shear rate or, equivalently, to small inelasticity. This allows us to go 
beyond the usual Navier-Stokes hydrodynamics. The hydrodynamic equations (|63fl - l|65|) are the starting point of the 
linear stability analysis of the USF of the next Section. 



V. LINEAR STABILITY ANALYSIS OF THE STEADY SHEAR FLOW STATE 



As said in the Introduction, computer simulations ; 32| have clearly shown that the USF state is unstable with 
respect to lon g enoug h wavelength perturbations. These results have been also confirmed by different analytical 
results j23, 134L ISfll l36j |. most of them based on the Navier-Stokes description that applies to first order in the shear 
rate. However, given that USF is inherently non-Newtonian 27], the full nonlinear dependence of the transport 
coefficients on the shear rate is required to perform a consistent linear stability analysis of the nonlinear hydrodynamic 
equations H63|) - (|65|) with respect to the USF state for small initial excitations. This analysis allows one to determine 
the hydrodynamic modes for states near USF as well the conditions for instabilities at long wavelengths. A growth of 
these modes signals the onset of instability, which is ultimately controlled by the dominance of nonlinear terms. Note 
also that while all the works have been mainly devoted to dense systems, much less attention has been paid to dilute 
gases. 

Let us assume that the deviations Sx^(r,t) — x^(r,t) — cco Al (r) are small, where fe M (r,i) denotes the deviation of 
{n, u, T} from their values in the USF state {no, uo, To}. The quantities in the USF verify 

Vn = VT = 0, u = a ■ r, d t T = 0. (66) 
Now, let us linearize Eqs. (|63|I - H65|I with respect to 

{fa M (r, t)} = {Sn(r, t), <5T(r, t), <5u(r, <)} . (67) 
The resulting set of five linearized hydrodynamic equations follows from Eqs. (|63[1 - H65|> : 

d 

d t 5n + ay—Sn + n ■ Su = 0, (68) 



3 d 

-n d t 5T + ay— ST + aS ix Su y + < 
2 Ox 



(d n pW)5n+(d T p£>)5T 



, r kl ~ ar lxykl 



Sn 
> 

'n 



-XonoT 2 



dSuk 
dri 
3 ST 
2To 



My 



dridrj 



K 



ry 

d 2 ST 
1 dri dri 



3 dSu k 



(69) 



d t Su k + ay—Su k 



aS kx Su v 



mno 



<f) p(°h dST 
-\°r p kt > — 



org oriTj 



d 2 S 



(70) 



Here, it is understood that the pressure tensor and its derivatives with respect to n and T, the cooling rate Co 
and the transport coefficients 77^7^, /iy , and are evaluated in the steady USF state. 

To analyze the linearized hydrodynamic equations (|68|) - (|70|l it is convenient to transform to the local Lagrangian 
frame, r[ — ri — ta^rj. The Lees-Edwards boundary conditions then become simple periodic boundary conditions in 
the variable r' 29] . A Fourier representation is defined as 



6x^,t)= / dr'e lk r fe„(r,t) 



J dr e^W-'&Epfoi), 



(71) 



where in the second equality ki(t) = kj(5ij — tdji). Periodicity conditions requires that ki = 2niTr/Li, where rij 
are integers and Li are the linear dimensions of the system. In this Fourier representation, the resulting set of five 
linear equations defines the hydrodynamic modes, i.e., linear response excitations to small perturbations. If at least 
one of the modes grows in time, the reference USF state is linearly unstable. Given the mathematical difficulties 
involved in the general problem, for the sake of simplicity, here I consider two kind of perturbations: (i) perturbations 
along the velocity gradient direction only (k x = k z = 0; k y 7^ 0) and (ii) perturbations in the vorticity direction only 
{k x — k y — 0; k z ^ 0). In both cases, the linearized hydrodynamic equations have time-independent coefficients. 
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A. Perturbations in the velocity gradient direction (k x — k z — 0; k y 7^ 0) 



Let us consider first perturbations along the y direction only. In this case, Eqs. (JBSJl-lJZDl hi this Fourier represen- 
tation can be written in the matrix form 



<9 T &r* + F^Sxl = 0, 
where the dimensionless quantities r = v^t and <5x* = {pk, @k, w^}, with 

Sn „ ST Su 



Pk 



) w K ry~\ ) *^k r— : ; 

n T y'To/m 

have been introduced. The matrix is 

F M „ = 2C8 ll2 5ui + C&^bvi + 0*^3^4 - ik*G,j, v + k* 2 H^ u , 



(72) 



(73) 



(74) 



where a* — a/vo, v$ is the collision frequency H21JI of the reference state and k* = i$k, £q = i/To/m/Vo being of the 
order of the mean free path. In addition, we have introduced the coefficient 



C(a) 



1 



a*(l + a*d:)P* 



and the square matrices 



(l-a*d a .)P* y {l~^a*d a ,)P* 

(i-a*d a ,)p; y (1 





1 a*d a .)R 



{■^xy ® Vxyxy) Ccy 3 i-^yy ® Vxyyy) ~^ Cyy 



y 
yy 













0/ 



(75) 



(76) 




















\ 




3^yy 


5 * 

3^yy 





















Vxyxy 


Vxyyy 















Vyyxy 


^yyyy 







V 














Vzyzy 


) 



have been also introduced. Here, P£ — /noT a and 



(77) 



C* =-^-a 2 )(Pkt-hl)v*kUy 



(78) 



The eigenvalues A AI (fc,a) of the matrix F(k,a) determine the time evolution of <5i*(fc,i). In the case that the real 
parts of the eigenvalues A M (fc, a) are positive, then the USF state will be linearly stable. Before considering the general 
case , it is convenient to consider some special limits. Thus, in the elastic limit (a = 1), the hydrodynamic modes 
of the Navier-Stokes equations (for the particular case considered here and in the context of the BGK model) are 
recovered ^l], namely, two sound modes, a heat mode and a two-fold degenerate shear mode. To second order in k* 
they are given by 



\^(k,a 




,*2 



- -k* 
3 



■„*2 



„*2 



(79) 



and consequently, excitations around equilibrium are damped. It is also quite illustrative to get the modes by setting 
k = 0, namely, consider small, homogenous deviations from the steady shear flow state. In this case, it is easy to see 
that pk and w Vt k are constant and 



w x .k{r) = w X}k (0) - a,TWy,k(0), 
0h(r) =6 k (0)e- Cr ~2 P k(0). 



(80) 
(81) 
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FIG. 5: Dependence of C(a) on the coefficient of restitution a. 




FIG. 6: Dispersion relations for a granular gas with a — 0.8 in the case of perturbations along the velocity gradient direction. 
Only the real parts of the eigenvalues is plotted. 

The mode associated with w x .k is unstable to an initial perturbation in w Vt k, leading to an unbounded linear change 
in time. However, stability is still possible at finite k if this behavior is modulated by exponential damping factors. 
With respect to the temperature field, initial disturbances decay at r — > oo if the coefficient C(a) > 0. Figure shows 
that the coefficient C is positive for any value of a and so, this mode is stable with a finite decay constant. 

The analysis for k ^ requires to get the eigenvalues A M (/c*, a) with the full nonlinear dependence of k*. However, 
the structure of F(fc, a) shows that the perturbation Sxt, oc 5u z is decoupled from the other four modes and hence can 
be obtained more easily. This is due to the choice of gradients along the y direction only. The eigenvalue associated 
with this mode is positive and is simply given by 

A 5 (fc, a) = v* zyzy k* 2 , rfzyzy = (ft + f»)2 ' ( 82 ) 

where C* is defined by Eq. J22). The remaining modes correspond to ph, Ok and the components of the velocity field 
w x .k and Wy^k- They are the solutions of a quartic equation with coefficients that depend on k* and a. The results 
show that Re A^(fc*,a) > for all the values of the coefficient of restitution a and consequently, the flow remains 
stable to this kind of perturbations. As an illustration, the dispersion relations for a gas with a = 0.8 are plotted in 
Fig. HO It is apparent that all the real parts of the eigenvalues are positive in the range of values of wavenumber 
k* considered. Our conclusion agrees with previous stability analysis (3^, based on the Navier-Stokes constitutive 
equations where it was found a minimum value of solid fraction (around 0.156) below which the USF is stable. Given 
that our system is a dilute gas (zero density), the present results confirm previous findings when one uses the improved 
transport coefficients. 

B. Perturbations in the vorticity direction (k x = k y = 0; k z ^ 0) 

The variation of the hydrodynamic modes with wavenumber k = k z in the vorticity direction is considered next. 
This situation has not been widely studied in the literature since most of the studies have been focussed on 2-D 
flows due to the relative computational efficiency with which they can be analyzed. Here, for the sake of simplicity, 
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FIG. 7: Stability lines k*(a) corresponding to the perturbation along the vorticity direction. The solid line corresponds to 
the results derived here while the dashed line refers to the results obtained from the Navier-Stokes approximation. The region 
above the curve corresponds to the stable domain, while the region below the curve corresponds to the unstable domain. 

I consider perturbations for which 8u x = Su y = and so, the eigenvalues X fJl (k*,a) obey a cubic equation. The 
analysis is similar to the one carried out in the previous section and so, details will be omitted. For a given value 
of a, it can be seen that this dispersion relation has one real root and a complex conjugate pair of damping modes. 
The instability arises from the real root since this mode Xu(k*,a) > if k* is larger than a certain threshold value 
k*(a). This value can be obtained by solving X ll (k*,a) = 0. As a consequence, the USF state is linearly stable 
against excitations with a wavenumber k* > k* (a) . The stability line k* (a) is plotted in Fig. [7] as a function of the 
coefficient of restitution. Above this line the modes are stable, while below this line they are unstable. For comparison, 
the corresponding stability line obtained from the approximations made in previous works |33l l35| is also plotted. 
This line can be formally obtained from the results derived in this paper when one replaces the expressions of the 
coefficients ijijkt, Kij, and fj,ij by their corresponding Navier-Stokes expressions It is apparent that the Navier- 
Stokes approximation captures the qualitative dependence of k* on a, although as expected quantitative discrepancies 
between both descriptions appear as the dissipation increases. Thus, for instance, for a — 0.8 the discrepancies 
between both approaches are about 22 % while for a = 0.5 the discrepancies are about 49%. The prediction of a 
long-wavelength instability for the USF state has been observed in early molecular dynamics simulations |32ll and 
qualitatively agrees with the previous analytical results based on the Navier-Stokes equations (3^, l34l l35l l3Gl|. At a 
quantitative level, the lack of numerical results from the Boltzmann equation prevent us to carry out a more detailed 
comparison to confirm the results derived from this kinetic model. We hope that the results offered here will stimulate 
the performance of such computer simulations. 

VI. SUMMARY AND DISCUSSION 

The objective of this paper has been to study the transport properties of a granular gas of inelastic hard spheres for 
the special nonequilibrium states near the uniform (simple) shear flow (USF). Although the derivation of the Navier- 
Stokes equations (with explicit expressions for the transport coefficients appearing in them) from a microscopic 
description has been widely worked out in the past 0, , the analysis of transport in a strongly shearing granular 
gas has received little attention due perhaps to its complexity and technical difficulties. Very recently, a generalized 
Chapman- Enskog method has been proposed to analyze transport around nonequilibrium states in granular gases . 
In the case of the USF state, due to the anisotropy induced in the system by the presence of shear flow, tensorial 
quantities are required to describe the momentum and heat fluxes instead of the usual Navier-Stokes transport 
coefficients 0, Q . In this paper we have been interested in a physical situation where weak spatial gradients of density, 
velocity and temperature coexist with a strong shear rate. Under these conditions, the corresponding generalized 
transport coefficients characterizing heat and momentum transport are nonlinear functions of both the (reduced) 
shear rate a* and the coefficient of restitution a. The determination of such transport coefficients has been the 
primary aim of this paper. 

Due to the difficulties embodied in this problem, a low-density gas described by the inelastic Boltzmann equation 
has been considered. Although the exact solution to the Boltzmann equation in the (steady) USF is not known, a good 
estimate of the relevant transport properties can be obtained by means of Grad's method [22, |23, |27j] . The reliability 
of this approximation has been recently assessed by comparison with Monte Carlo simulations of the Boltzmann 
equation [22L l2Sj . Assuming that the USF state is slightly perturbed, the Boltzmann equation has been solved by 
a Chapman-Enskog-like expansion where the shear flow state is used as the reference state rather than the local 
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equilibrium or the (local) homogeneous cooling state. Due to the spatial dependence of the zeroth-order distribution 
j(o) (reference state), this distribution is not in general stationary and only in very special conditions has a simple 
relation with the (steady) USF distribution 0. Here, since one the main goals has been to address a stability analysis 
of the USF state, for practical purposes my results have been specialized to the steady state, namely, when the 
hydrodynamic variables satisfy the balance condition l|52[) . In this situation, the (reduced) shear rate a* is coupled 
with the coefficient of restitution a [see Eq. I(25|) ] so that the latter is the relevant parameter of the problem. In the 
first order of the expansion the momentum and heat fluxes are given by Eqs. I|47|) and l|48|) . respectively, where the 
set of generalized transport coefficients rjijki, f^ij, and Ky are given in terms of the solutions of the linear integral 
equations lf53jl - lfST|) . As expected, there are many new transport coefficients in comparison to the case of states near 
equilibrium or cooling state. These coefficients provide all the information on the physical mechanisms involved in 
the transport of momentum and energy under shear flow. 

Practical applications require to solve the integral equations l|55|) - (|57|l . which is in general quite a complex problem. 
In addition, the fourth-degree velocity moments of USF (whose evaluation would require to consider higher-order terms 
in Grad's solution 1)18(1 of the Boltzmann equation) are needed to determine the coefficients Kij and jMj- To overcome 
such mathematical difficulties, here a kinetic model of the Boltzmann equation [3lT | has been used. This kinetic model 
can be considered as an extension of the well-known BGK equation to inelastic gases. Although the kinetic model is 
only a crude representation of the Boltzmann equation, it does preserve the most important features for transport, 
such as the homogeneous cooling state and the macroscopic conservation laws. The model has a free parameter (3 to 
be adjusted to fit a given property of the Boltzmann equation. Here, [3 is given by Eq. (|23[) to get good quantitative 
agreement of the Navier-Stokes shear viscosity coefficient obtained from the Boltzmann equation. Furthermore, this 
choice yields the same results for rheological properties in the USF problem as those derived from the Boltzmann 
equation by means of Grad's method. On the other hand, given that the model does not intend to mimic the behavior 
of the true distribution function beyond the thermal velocity region, discrepancies between the kinetic model and the 
Boltzmann equation are expected beyond the second-degree velocity moments (which quantify the elements of the 
pressure tensor). Nevertheless, a recent comparison with Monte Carlo simulations of the Boltzmann equation [2^] 
have shown the accuracy of the kinetic model predictions for the fourth-degree moments. As illustrated in Fig.^ the 
semi-quantitative agreement between theory and simulation is not restricted to the quasielastic limit (a ~ 0.99) since 
it covers values of large dissipation (a > 0.5). The use of this kinetic model allows one to get the explicit dependence 
of the generalized transport coefficients on the coefficient of restitution. This dependence has been illustrated in some 
cases showing that in general the deviation of the transport coefficients from their corresponding elastic values is quite 
significant. 

With these new expressions for the fluxes, a closed set of generalized hydrodynamic equations for states close to USF 
has been derived. A stability analysis of these linearized hydrodynamic equations with respect to the USF state have 
been also carried out to identify the conditions for stability in terms of dissipation. Two different kind of perturbations 
to the USF state has been analyzed: (i) perturbations along the velocity gradient only (fc„ ^ 0) and (ii) perturbations 
along the vorticity direction only (k z ^ 0). In the first case, previous results 133. |35l have shown that the USF is 
stable for a dilute gas while the USF becomes unstable in the second case for all a [42] • These results agree with these 
findings and the USF is unstable for any finite value of dissipation at sufficiently long wave lengths when disturbances 
are generated in the orthogonal direction to the shear flow plane. On the other hand, as expected, quantitative 
discrepancies between our results and those given |33l l35| from the Navier-Stokes approximation become significant 
as the dissipation increases. These differences have been illustrated in Fig. [7| for the stability line. Although the 
instability of the USF has been extensively studied for many authors by using a Navier-Stokes description |33ll33 .l35| 
as well as solutions of the Boltzmann equation in the quasielastic limit pa . l37j . I am not aware of any previous 
solution of the hydrodynamic equations where the generalized transport coefficients describing transport around USF 
were taken into account. The analytical results found in this paper allows a quantitative comparison with numerical 
solutions to the Boltzmann equation for finite dissipation. As happens for the USF problem for elastic [2^, 0. 
and inelastic I^H^ gases, one expects that the results reported here compare well with such simulations, confirming 
again the reliability of the kinetic theory results to characterize the onset and the first stages of evolution of the 
clustering instability. We hope to carry out these simulations in the next future. 

On the other hand, the stability analysis performed here has only considered spatial variations along the y and z 
directions. More complex dynamics is expected in the general case of arbitrary direction for the spatial perturbation. 
This will be worked elsewhere along with comparison with direct Monte Carlo computer simulations of the Boltzmann 
equation. Another possible direction of study is the extension of the present approach to other physically interesting 
reference states, such as the nonlinear Couette flow. This is a more realistic shearing problem than the USF state 
since combined heat and momentum transport appears in the system. Given that an exact solution to the kinetic 
model used here is known for the Couette flow problem |3(j, the reference distribution for the Chapman-Enskog-like 
expansion is available. 
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APPENDIX A: CHAPMAN-ENSKOG EXPANSION 

Inserting the expansions l|32|l and (|34|1 into Eq. Ij26(l . one gets the kinetic equation for 



d t (0) -aV v ^-+C) f& = a?* + (V + u ) • V /W, (Al) 







v dV x 

where C is the linearized Boltzmann collision operator 



£X=-(j[/(°),X] + J[X,/(°)]). (A2) 

The velocity dependence on the right side of Eq. IjAlfl can be obtained from the macroscopic balance equations to 
first order in the gradients. They are given by 

d[ 1] n + u • Vn = -V • (nSn), (A3) 



5 t (1) <5u+ (u + <5u) • X75u = — V ■ P (0) , (A4) 

P 

^a t (1) T+^n(uo + <5u)-Vr + aP« + p(°) : V<5u = -\pt (1) , (A5) 
where p = mn is the mass density, 

P[p = f dcmc iCj f^(c), (A6) 



and 



C (1) = — / dcmc 2 Cf (1) . (A7) 
3pJ 



Use of Eqs. l|A"3 )l -l|A"5 |l in Eq. (JJTJ yields 



df ] -aV y 4r+ - C {1) T^- = Y„ • Vn + Y T ■ VT + Y u : V<Ju, (A8) 



where 



dn 1 p dSuj dn 



df (o) ld f(0) dP (® 
YT,i = -^rCi + ~— (A10) 



dT 1 p dSuj dT ' 

5/(0) 2 



07^ oj^ a o/^ / (Q) \ 
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According to Eqs. ljA9(l - l|A10|) . Y u ,ij has the same symmetry properties (|16|) as the distribution function /(°) while 
Y ni and Yt,i are odd functions in the velocity c. 
The solution to Eq. I|A8(I has the form 



= X n>i (c)Vin + X T ,i(c)ViT + \„,,ic)V,du l 



(A12) 



Note that in Eq. (jAl If) the coefficients r\ijki are defined through Eq. I|49|) . Substitution of the solution (| Al 2|) into the 
relation l|A7(l allows one to write the cooling rate in the form 



C (1) = Cn,Nin + Cr.iViT + (ujiViSu 



'3 ' 



(A13) 



where 



Cn,z 



Hp 



dcmc 2 C I Xt } 



(A14) 



X 



However, given that A^ and Aj^i are odd functions in c [see for instance, Eqs. (|A19(1 and i|A20|l below], the terms 
proportional to Vn and VT vanish by symmetry, i.e., 



Cn,i — Cr,i — 0. 

Thus, the only nonzero contribution to Q 1 ' comes from the term proportional to the tensor ViSuji 

C (1) = Cu jiViSuj. 



(A15) 



(A16) 



An estimate of the tensor Q u ,ij has been made in Appendix[5]by considering the leading terms in a Sonine polynomial 
expansion of the distribution / W. Its expression is given by Eq. I|B9|1 . As expected, ( u ij vanishes in the elastic limit 
(« = !)• 

The coefficients X n ^, Xta, and X u ^j are functions of the peculiar velocity c and the hydrodynamic fields. In 
addition, there are contributions from the time derivative <9 t acting on the temperature and velocity gradients given 
by 



sfViT = v,a t (0) T 



* (1 - „ a „,P iS ) - v,„ - (Iw© + 1<<». ) v,t, 



(A17) 



(0), 



jkViOUk- 



Substituting (|A16|I into l|A8|l and identifying coefficients of independent gradients gives the set of equations 



T 

X n a H — 
n 



2a 



g(l-m)„)Pj«>-C (0) 



At,i — 



(A18) 



(A19) 



J^i? + ^C (0) ) a T + |r(^) + §<<°> + ac v A - £ 



At,i — 



(A20) 



2 

3ri 



A 



u,M — aoky 



Upon writing Eqs. (|A19|) - (|A21|I . use has been made of the property 



d^X = 



dX 



a (0) T 



dT * 
2 



dX 



85- 1 



(A21) 



(A22) 



where in the last step we have taken into account that X depends on Su through c = V <5u. 
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APPENDIX B: EVALUATION OF THE COOLING RATE 

In this Appendix the contribution C,u,ij to the cooling rate is evaluated by expanding X Ut ij as series in Sonine 
polynomials and taking the lowest order truncation. The tensor C, U} ij is given by 

Cu,ij = ^- J dci mc\LX u sj 

= ~J dc 1 mc 2 {j[c 1 \f ( -°\X u>ij ] + J[ Cl \X Utij ,fW]}. (Bl) 

A useful identity for an arbitrary function h(ci) is 

J dc 1 h(c 1 )J[c 1 \f,g]=a 2 J dci J dc 2 f( Cl )g(c 2 ) J da Q(a ■ g)(ff • g) [h(c'{) - h(a)] , (B2) 

where g = Ci — c 2 and 

ci' = ci-i(l+a)(£-g)3. (B3) 



Using l|B2|) . Eq. (JBlfl can be written as 



Cu,ij = |V(1 - a 2 ) J dc x J dc 2 f {Q) ( Cl )X u ^(c 2 ) J d*Q{S ■ g)(S • g) 3 . (B4) 



The integration over a in i|B4|) yields 



77? 



(u,ij = y^^ 2 (1 " a 2 ) J dc x J dc 2 / (0) ( Cl )A M , y (c 2 ). (B5) 

This equation is still exact. To perform the integrals over ci and c 2 one takes the Grad approximation (|18|l to /W 
and expands -X" Ui jj in Sonine polynomials. In this case and according to the anisotropy of the USF problem, one takes 
the approximation 

X aM (c) -»■ -■^j^DijrjijkiMc), (B6) 



where 



is the Maxwellian distribution and 



A(e) .„(^^(_-) (B7) 



1 



Aj(c) = m I CjCj - -c <%J . (B8) 

Next, change variables to the (dimensionless) relative velocity g* = (ci — c 2 )/vo and center of mass G* = (ci+C2)/2tjo, 
where vq = y / 2T/m is the thermal velocity. A lengthy calculation leads to 

Cu,ij = "^(l-« 2 ) / dg* J dG* g^e- 2G ' 2 e^' 2 / 2 



G* k G* t G* m G* n - —g* 2 G* 2 {6 km S en + S kn S em ) + —g k gigm9*, 



P, 



in II 



nT 



1 a 2 ,/^(l - a 2 ) ( % - * w ) (B9) 



15 V toT V nT 

Of course, when a = 1, then ( u jj = 
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APPENDIX C: BEHAVIOR OF THE ZEROTH-ORDER VELOCITY MOMENTS NEAR THE STEADY 

STATE 



This Appendix addresses the behavior of the velocity moments of the zeroth-order distribution / W near the steady 
state. Let us start with the elements of the pressure tensor pfj . In the context of the Boltzmann equation and by 
using Grad's approximation (|18J) . they verify the equation 

~ {l aP *y + rc(0) ) ^ + + = v V ~ pS -) + CP ^} • (C1) 

Since we are interested in the hydrodynamic solution, the temperature derivative term can be written as 

Td T Pf = Td TP P* = p (l _ i " A W , (C2) 



v >■* v * y 2 da* J V} 

where P*j — Py /p. Upon deriving l|C2(l . use has been made of the fact that the dimensionless pressure tensor P.* 
depends on T only through its dependence on the reduced shear rate a* = a/v(n,T). In dimensionless form, the set 
of equations (|C1|) become 

- Qa*P* y + C*) (l - l a *^) p *i + a te p J* + a h p u = -[13 (Ai - Sy) + CP*} , (C3) 



where 



t(0) 5 

C = i- = -(l- a 2 ). (C4) 



Let us consider the elements P* y and P* y = P* z . From Eq. (|C1|I . one gets 

\^K V + C*) (l - K V + a*PZ y = - 08 + O Ply, (C5) 

- (|a*i^ + C*) (l - \a*-^j P; y = -(/? + C)P; y + (3- (C6) 

This set of equations have a singular point corresponding to the steady state solution, i.e., when a*(T) — a* where 
a* (a) is the steady state value of a* given by Eq. I|25|l . Since we are interested in the solution of Eqs. I|C5|1 and (|C6|) 
near the steady state, we assume that in this region P* y and P* y behave as 

(dP* \ 



[dP 



where the subscript s means that the quantities are evaluated in the steady state. Substitution of 1|C7(I and (|C8|) into 
Eqs. ljC5|) and i|C6(l allows one to determine the corresponding derivatives. The result is 

r)p* \ n *C 4- P* 

ur vv_\ - ap* s x y> s fro} 



da* J s vv ' s 2af C + 6/3 + 3C 
where C = (dP* y /da*) is the real root of the cubic equation 

2a* 4 C 3 + 12af(C + f3)C 2 + ^(7C* 2 + U(*/3 + 4/3 2 )C + 9/3(C* + /3)" 2 (2/3 2 - 2C* 2 - PC)- (CIO) 

Equations (|C9() and IjClOf) can be also obtained from a different way. Let us write the set of equations (|C5|I and (|C6|) 

as 

f)P* —IP* —P* (B — -P*a*) 

Ur xy _ zr yy a* r xy V J Z r xy a j fCUl 

8a* C + h*P* 
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®Pyv _ ^yy \P 3^v 



da* ' a* (C* + | a*P* 



(C12) 



xyj 

In the steady state limit (a* — > a*), the numerators and denominators of Eqs. I|C11|I and l|C12(l vanish. Evaluating 
the corresponding limit by means of l'Hopital's rule, one reobtains the above results l|C9|l and l|C10|) . This procedure 
can be used to get the behavior of the remaining velocity moments near the steady state. 

The behavior of the fourth-degree velocity moments of the distribution /(°) near the steady state is also needed to 
determine the transport coefficients fiij and associated with the heat flux in the first-order solution. To evaluate 
this behavior we use the Boltzmann kinetic model (|59() . Let us introduce the velocity moments of the zeroth-order 
distribution 

Kmm= J *<£cfcc*"/ (0) (c) (C13) 

These moments verify the equation 



^i,/..'-*y^*' |C14) 



where k = k± + k 2 + k^ and N klt k 2 ,k s are the velocity moments of the Gaussian distribution fg. As before, the 
derivative 9rM^ fe fca can be written as 



2T\ k/2 1 



n 

m J 



(k-a*d a *)M* kuk2M (a*). (C15) 



In dimensionless form, Eq. I|C14|I become 

- (l a * P *y + C*) \( k - ^d a ,)M* kiMM + k ia *M* ki _ 1M+lM + (/} + \OM* kl 



k2,k3 

«,fe 2 ,%= 5 ( C16 ) 



where N£ k2 k3 are the reduced moments of the the Gaussian distribution given by 

'ki + 1\ fk 2 + 1 



NL ,,,,, = -- 3/2 r ^ ^ ^ (Ci7) 



if fci, k 2 , and k% are even, being zero otherwise. Equation (|G16|I gives the expressions of the reduced moments 
M% k k a in the steady state. To get d a *M ki k2 k3 in the steady state, we differentiate with respect to a* both sides 
of Eq. 1C16|) and then takes the limit a — > a*. In general, it is easy to see that the problem becomes linear so that 
it can be easily solved. To illustrate the procedure, let us consider for simplicity the moment Mq 40 , which obeys the 
equation 



\**P; y + Cj (2 - \**da^j Af * 40 + 09 + 2C) M* 40 -\p = 0. 



(C18) 



From this equation, one gets the identity 



2 - \a*d a , ) M 



[P* xy + a*{d a ,P* xy )} [2-\a*d a ,)M { 



041) 



3 

+(/? + 2O9 o .M * 40 = (C19) 

In the steady state limit, Eq. I|52|l applies and the first term on the left hand side vanishes. In this case, one easily 
gets 

fa' M ™J S - a* sXs +20 + <* 



^ M 040 = -r , Z i ^ M oV s , (C20) 
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where 

2 

is a known function and 



P*. . + a 



xy,s s \ Qa 



dP* 



(C21) 



Proceeding in a similar way, all the derivatives of the form d a *M* can be analytically computed in the steady state. 

APPENDIX D: KINETIC MODEL RESULTS IN THE STEADY STATE 

In this Appendix, I display the results obtained from the model kinetic equation chosen here for the determination 
of the generalized transport coefficients. In the model, the Boltzmann collision operator is replaced by the term [3l| 

J[/,/]->-M/-/o) + ~-(c/), (Dl) 

where v and (3 are given by Eqs. i|21|) and l)23|) . respectively, /o is the local equilibrium distribution (|19|l and ( is the 
cooling rate fl2"|) . 

1. Steady state solution for the (unperturbed) USF 

Let us consider first the steady state solution to the (unperturbed) USF problem. In this case, <5u = and so 
c = V. The one-particle distribution function /(V) obeys the kinetic equation 

-aV v ^r/(V) = -fivtf -fo + ^-^L. (V/) , (D2) 

where here £ has been approximated by its local equilibrium approximation given by Eq. I|22|) . The main advantage 
of using a kinetic model instead of the Boltzmann equation is that the model lends itself to an exact solution 0, l28j . 
It can be written as 

/(V)=n(^) ^ /*(<■), ^ = ^V, (D3) 

where the reduced velocity distribution function /* is a function of the coefficient of restitution a and the reduced 
peculiar velocity 

r°° — r — 

/*(£) = it- 3 / 2 dae-^-tOexp -e Cs (£ + sa • £) 2 . (D4) 



Here, a = z/{y0) and C = C^°V( jy /3)- ^ n as been recently shown that the distribution function l|D4|l presents an 
excellent agreement with Monte Carlo simulations in the region of thermal velocities, even for strong dissipation [28| . 

The explicitly knowledge of the velocity distribution function allows one to compute all the velocity moments. We 
introduce the moments 

M kM * = J dvV x k >V y k *V z k >f(V) (D5) 

According to the symmetry of the USF distribution (|D4fl . the only nonvanishing moments correspond to even values 
of k\ + hi and k^. In this case, after some algebra, one gets [28| 



fc/2 

m I 



'2T\ ' 

M kuk2ik3 =n{—) M* kMa , (D6) 



23 



where the reduced moments k2 ks are given by 



7 \ ; i 



g = 

q-\-k\ —even 



(fci - 9)! 



r ' fci-g+l \ r / fc 2 +g + l \ r /fe 3 + 1 



(D7) 



with a — a/(vf3) — a*//3. It is easy to see that the expressions for the second degree-degree velocity moments 
(rhcological properties) coincide with those given from the Boltzmann equation by using Grad's approximation, Eqs. 



2. Transport coefficients 



Let us now evaluate the generalized transport coefficients rjijke,Kij, and fjtij in the steady state. They can be 
obtained from Eqs. I)55(l -I(57 () with the replacement given by Eqs. 1)61(1 and lfi2|l . With these changes, Eqs. 1)55 (1 - 1)57(1 
become 



dc x 



3 n 



(D8) 



9 1 \ 

Cy ~dc~ ~ 3° ( P * y ~ a * da * p xy) + C ) Xt ' 1 = Yt ' 1 ' 



(D9) 



d C (0) d \ 1 



dc (J )+ dT 1 



- aSj y X UtX e = Yuje. (D10) 
In order to get the transport coefficients Kij , fXij , and rjijkt , it is convenient to introduce the velocity moments 

(Dll) 



A (i) 



dcc kl c k2 c k3 X ■ 

x y z 



B 



ki,k 2 ,k 3 



J dct£c*'(£x T ,i, 



(D12) 



L, fe 1 ,fc 2 ,fc 3 



tjc r kl r k2 r k3 X 



(D13) 



The knowledge of these moments allows one to get all the transport coefficients of the perturbed USF problem. Now, 
we multiply Eqs. ((D8(1 - I(D10() by c kl c k l 2 c k3 and integrate over velocity. The result is 



(D14) 



(D15) 



au 3V^ki,k 2 M — / ULC a; H, c z '"J 



(D16) 
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Here, M ( k 0) k k are the moments of the zeroth-order distribution and we have introduced the quantities 



The right-hand side terms of Eqs. l|D14|l - l|D16Jl can be easily evaluated with the result 



(D17) 



A 



k±,k 2 ,k 3 



dcCg 1 C y 2 C z 



idP, 



(<>) 



x (5j x kiM kl _i tk2ik3 + 5j y k 2 M kuk2 _i jk3 + 5j Z k 3 M klMM -i) 



2TY 
m J 



(1 - a*d a ,)M^ i+StxM+StvM+Slx - ~(1 - a*d a *)Pl 



x (S^hM*^^ +S jy k 2 M* iM _ 1M +^A: 3 M fc * lifc2ifc3 _ 1 )] , 



(D18) 



fel,fe2,fe3 



^j/ ~z 



<9 



1< } 



Qr^Mki+Sl x ,k2+Sly,k3 + &ez ^ Qrp 

x (^ K fciM fcl _i :fc2 ^ 3 + 5jyk2M kuk2 -i, k3 + 8 Jz k 3 M klM ^ 3 _i) 



2T 



^(k + 1 - a*9 a .)Af fe * 1+(5fxifc2+5f!(ife3+( 5 fz - ^(1 - ^a*d a ,)P;, 



x fexfciM^^^^ +^feM fc * i fc2 „ l fc3 +S j zk 3 M k ¥ i k2 k3 _ 1 )] 



(D19) 



C 



k 1 ,k 2 ,k 3 



J dcc^c^c k /Y u , je 



1 — n- 



d_ 

dn , 



dT 



M, 



k 1 ,k 2 ,k 3 



-M kljk2tka (S jx Si> x ki + 5 jy 5iyk 2 + 5 JZ 8i z k 3 ) 
—ki8j X (5gyM kl -i tk2+ i tk3 + (5£ Z Mfe 1 _i : fc 2: fc 3+ i) 
—k 2 Sj V [^ix^lki+x.ki-^M + Se z Mk ly k 2 -i,k 3 +i) 
—k 3 Sj Z (5g x M kl+ i^ k2 . k3 -i + SiyM^^+iM-i) 

f2T\ h/2 
~ n \m) i- S ^ a * d ** M Lk 2 ,k 3 

~3^T ~ aVx y^) ^ ~ a * d ^ M kiM,k 3 

+ M k 1: k 2 ,k 3 (SjxStxki + 5 ]y 5iyk 2 + S JZ 8t z k 3 ) 
+k 1 S jx (S ey M^_ ltk2+1>k:3 + S ez M^_ 1MM+1 ) 
+k 2 S jy (S ix Ml 1+lM _ hk3 + S ez M k * iM _ 1M+1 ) 
+k 3 S jz (5 tx Ml i+XMM _ x + h y Ml xM+1M _ x ) 



(D20) 



Here, M ki k2 k are the reduced moments of the distribution /W defined by Eq. I|D6|) . In the steady state, M£ fc2 k is 
given by Eq. I|D7|I while the derivatives d a *M ki k2 k3 can be obtained by following the procedure described in Appendix 



The solution to Eqs. pTU|) - (lD16|) can be written as 



A 



(i) 

ki,k 2 ,k 3 



fci 



9=0 



fcC 

2 



-(1+9) 



A'i! 



(fci-g)! 



■^kx-qM+qM WnD k x -q,k 2 +lM 



,(>) 



(D21) 
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B 



(i) 

ki,k 2 ,k 3 



■7^ -(1+9) 



■ W T + — 



9=0 



(fcl " q) 



I ^fei -q,k 2 +q,k 3 ' 



(D22) 



ki,k 2 ,k 3 



fcl 



9=0 



fc_C 

2 



T\ -(l+<?) 



(fci - 9)! 



. ^df) 1 / 2T\ 



'k 1 —q,k 2 +q.k 3 ""JJ/^fci — g,/S2+g,fc3 2 



fc/2 



TO / 



(D23) 



where = ujt/{v[3). From Eqs. I)D21|) - (|D23|) one can get the expressions for the transport coefficients fey, /xy, and 
T7ijfc£ in terms of /?, £ and a. 
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